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Abstract 

By using a semi-analytical dynamical mean-field approximation previously pro- 
posed by the author [H. Hasegawa, Phys. Rev. E , 70, 066107 (2004)], we have 
studied the synchronization of stochastic, small-world (SW) networks of FitzHugh- 
Nagumo neurons with diffusive couplings. The difference and similarity between 
results for diffusive and sigmoid couplings have been discussed. It has been shown 
that with introducing the weak heterogeneity to regular networks, the synchro- 
nization may be slightly increased for diffusive couplings, while it is decreased for 
sigmoid couplings. This increase in the synchronization for diffusive couplings is 
shown to be due to their local, negative feedback contributions, but not due to the 
shorten average distance in SW networks. Synchronization of SW networks depends 
not only on their structure but also on the type of couplings. 
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1 INTRODUCTION 



In recent years, much attention has been paid to complex networks such as small-world 
(SW) and scale-free (SF) networks (for a review see [1, 2]). The SW network, which has 
been proposed by Wattz and Strogatz (WS) [3, 4], is characterized by a large clustering 
coefficient and a small average distance. The original WS-SW network has been created by 
introducing the finite-degree heterogeneity to the regular network by random rewirings of 
links. Newman and Wattz (NW) [5] have proposed an alternative SW network, randomly 
adding shortcut links to the regular networks without rewirings. In SW networks, the 
degree distribution P(k) for a node to have k coupled neighbors has an exponential tail 
for a large k. In contrast, the degree distribution in SF networks, which was first proposed 
by Barabasi and Albert [6], is given by the power law as P{k) ~ A;~ 7 with the index 7 
(= 2 ~ 4). Since Barabasi and Albert [6] have proposed a growing SF network with 
preferential attachments of nodes, many models and mechanisms have been proposed 
not only for growing but also for non-growing SF networks with geographical and non- 
geographical structures. 

The interplay between structure and dynamics has attracted a great deal of atten- 
tion to the synchronization in complex networks. The synchronization in SW networks 
consisting of spiking neurons has been studied [7]- [15] with the use of Hodgkin-Huxley 
(HH) [7] [8], FitzHugh-Nagumo (FN) [9][10] Hindmarsh-Rose (HR) [11], integrate-and-fire 
(IF) [8] [12], and phase models [13, 14]. By using a more general class of models, dynam- 
ical properties including synchrony in SW and SF networks have been also investigated 
[16, 17, 18]. It has been, however, controversial whether the synchronization in complex 
networks is better or worse than that in regular networks. Most of calculations have 
shown that the synchronization in SW networks is better than regular networks because 
of the shorten average distance in the former [7, 11, 13, 14, 15, 16]. On the contrary, it has 
been shown that the average distance is not necessarily correlated with the synchroniz- 
ability of the networks [10] [17]. Some have claimed that the synchronization is increased 
or decreased depending on the adopted parameters or calculation conditions [8] [9] [12]. 

In a previous paper [10], we have developed a semi-analytical theory for SW networks, 
by generalizing the dynamical mean-field approximation (DMA) which was originally pro- 
posed for regular networks with all-to-all couplings [19, 20]. The method newly developed 
in [10] is applicable to SW networks with a wide range of couplings covering from local 
to global and/or from regular to random ones. In [10], we have taken into account three 
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kinds of spatial correlations: on-site correlation, the correlation for a coupled pair, and 
that for a pair without direct couplings. Our method has been applied to SW FN neural 
networks with sigmoid couplings [10], in which a coupling to the neuron i is given by 

lt\t) = J^G(x ] (t)), (1) 
j 

where J denotes the coupling strength, G(x) = 1/[1 + exp(— (x — 0)/a)\ is the sigmoid 
function with threshold 9 and width a, and the adjacent matrix is = Cji = 1 for 
a coupled (i, j) pair and otherwise. Calculations by DMA and direct simulations have 
shown that when random links are added to regular networks, the synchronization is 
decreased because of the introduced heterogeneity in SW networks [10]. 
Besides the sigmoid coupling of Eq. (1), the diffusive coupling given by 

lt\t) = Kj2c lJ H(x 3 (t)-x t (t)), (2) 

3 

has been widely employed for theoretical study on neural networks, where K stands for 
the coupling strength and H(x — y) the coupling function. Equations (1) and (2) model 
chemical and electrical synapses, respectively. Both chemical and electrical synapses 
exist in neocortex. Chemical synapses use a chemical neurotransmitter that is packaged 
presynaptically into vesicle, released in quantized amount, and binds to postsynaptic 
receptors. In contrast, electrical synapses are simpler in structure and function. They 
provide a direct pathway that allows ionic current to flow from the cytoplasm of one cell to 
that of another [21]. Although chemical synapses are by far the most abundant, electrical 
synapses also play an important role in neocortex. The purpose of the present paper is 
to apply the DMA to SW neural networks of FN neurons with diffusive couplings, and to 
compare the results for diffusive couplings to those for sigmoid couplings in [10]. This is 
expected to provide some insight to unsettled issue on the effect of the heterogeneity on 
the synchronization in SW networks mentioned above. 

The paper is organized as follows. In Sec. II, we have derived differential equations 
(DEs), applying the DMA to SW FN networks with diffusive couplings in order to trans- 
form the original stochastic DEs to deterministic DEs. In Sec. 3.1 and 3.2, we have 
reported numerical calculations for regular and SW networks, respectively. The final Sec. 
IV is devoted to conclusion and discussion. 
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2 Small-world networks of FN neurons 



2.1 Adopted model and method 

We have assumed that iV-unit FN neurons are distributed on a ring with the average 
coordination number Z and the coupling randomness p. Dynamics of a single neuron % in 
a given SW network is described by the non-linear DEs given by 

dx H (t) 



dt 
dx 2 j{t) 

dt 



with 



= F[x u {t)}-cx 2l {t) + l\ c) (t) + l\ e) {t)+Ut), (3) 
= b xu(t) - d X2i(t) + e, (i = 1 to N) (4) 

lt\t) = tf£ctf *!<(*)), (5) 

3 

l\ e \t) = AQ{t-t m )Q{t m + t w -t). (6) 

In Eq. (3)-(6), F[x(t)\ = k x(t) [x(t) - a] [1 - x(t)\, k = 0.5, a = 0.1, b = 0.015, d = 0.003 
and e = [19] [22] [23]: x Vi and rr 2 i denote the fast (voltage) variable and slow (recovery) 
variable, respectively: H(x) stands for the diffusive-type coupling: the adjacent matrix 
given by Qj = Cji = 1 for a coupled (i, j) pair and zero otherwise, self-coupling terms 
being excluded (cu = 0). By changing Z value, our model given by Eqs. (3)-(6) covers 
from local couplings (Z <^ N) to global couplings (Z — N — 1). We should, however, 
keep in mind that the electrical synapses by nature can only be produced among close 
neurons. The response of neuron networks has been studied to an external, single spike 
input given by lj e \t) with magnitude A and spike width t w applied for t in < t < t in + t w , 
Q(x) being the Heaviside function. Added white noises ^(t) are given by 

<Ut)> = 0, (7) 
<Ut)Cj(t')> = P 2 S tJ 5(t-t'), (8) 

where the average of < U(z,t) > for an arbitrary function of U(z,t) is given by 

<U(z,t)>= J.. .J dzU(z,t) Pr(z), (9) 

Pr(z) denoting a probability distribution function for 2iV-dimensional random variables 
z = ({x Ki }). 
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Our WS-SW network has been made after [3]. Starting from a regular network, N c h 
couplings among NZ/2 couplings are randomly modified such that Cj 3 - = 1 is changed to 
Cij = or vice versa. The coupling randomness p is given by 

2N ch 

p= Jtz> (10) 

which is and 1 for completely regular and random couplings, respectively. 

In DMA [19] [10], we will obtain equations of motions for means, variances and co- 
variances of state variables. Variables spatially averaged over the ensemble are defined 
by 

X M = ^ E K = 1 > 2 (11) 

i 

and their means by 

lM K (t) = ((X K (t))) c , (12) 

where the bracket < • > c denotes the average over the coupling configuration. As for 
variances and covariances of state variables, we consider three kinds of spatial correlations: 
(1) on-site correlation (7), (2) the correlation for a coupled pair ((), and (3) that for a 
pair without direct couplings (77): 

f 7 KjA , for % = j 
((5x Ki 8x Xj )) c = < C«,a, for i ^ j, Cij = 1 (13) 
[ r] KjX , for % ^ j, dj = 

where k, A = 1, 2 and 

5x Ki {t) = x Ki (t) - fi K {t). (14) 
In Eq. (13), j KiX , ( K ,\ and r] KyX are defined by 

7«,a(0 = (^E(^W^(t))\ , (15) 



N 

Wz 



CAt) = (^EE^Wl^f). (i6) 



^,A(t) = ^ Ar(Ar _ 1 z _ 1) E E(l " % - %) 5x Xj (t))^ . (17) 

For a later purpose, we define also the spatially-averaged correlation given by 

p K ,x{t) = I^EEM)^))). (is) 

= «<5X K (t)5X A (t)» c , (19) 
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where SX K (t) = X K {t) — (J, K (t). It is noted that 7 Kj a, Ck,a an d p Kt \ are not independent, 
obeying the sum rule given by 

Np K , x = 7k ,a + Z( K , X + (N-Z- l) VKtX . (20) 

In order to derive Eqs. (15)-(20), we have employed the decomposition: 

i = Sij + ii -<y[cij + (i -Cij)}, 

= $ij + Cij + (1 - Sij - 0^), (21) 

with ca — 0. 

In calculating means, variances and covariances given by Eqs. (12) and (15)-(18), 
we have assumed that (1) the noise intensity is weak, and (2) the distribution of state 
variables takes the Gaussian form. By using the first assumption, we expand DEs given 
by Eqs. (3)- (6) in a power series of fluctuations around means. The second assumption 
may be justified by some numerical calculations for stochastic FN [23] [24] and HH neuron 
models [25] [26]. It has been shown that for weak noises, the distribution of a membrane 
potential of a single FN or HH neuron nearly follows the Gaussian distribution, although 
for strong noises, the distribution deviates from the Gaussian, taking a bimodal form. 

Before closing Sec. 2.1, we briefly summarize the introduced variables and their mean- 
ings as follows: N, the number of neurons: Z, the average coordination number: p, the 
coupling randomness: K, the coupling strength: c^, the adjacent matrix: X K , the spa- 
tially average of the fast (k = 1) and slow (k = 2) variables; /i K , a mean value of X K ; 
7k,a, Ck,a, and i] K ^, the correlations of on-site, a coupled pair, and an uncoupled pair, 
respectively. 

2.2 Equations of motions 

We will obtain equations of motions for fi K (t), 7«,a(£), Ck,a(0> Vn,\(t) an d Pk,a(0- Readers 
who are not interested in mathematical details, may skip to Sec. 2.3, where our theoretical 
results are summarized. 

After some manipulations, we get the following DEs (the argument t being suppressed; 
for details, see appendix A): 

^ = fo + /271.1 - Cfi2 + Iext, (22) 

^ = b^-d^ + e, (23) 
dt 
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with 



dr): 



2,2 



dt 
drji j2 
dt 



2(a 7 i,i - c 7 i, 2 ) + 2KZh 1 (( 1 ,i - 7 i,i) + (24) 
2(^1,2-^2,2), (25) 
& 7 i,i + (a - d) 7 i j2 - c 72 ,2 + K^i(Ci,2 - 7 i,i), (26) 



#1,1 
dt 

d 7 2,2 

dt 
dt 

^ = 2(a Plil - c Pl)2 ) + ^, (27) 

dp 2 , 2 

dt 

dt 
<*Ci,i 
dt 

#2,2 

(it 

#1,2 
dt 

<fyl,l 

dt 



= 2(6pi i2 -dp2 i2 ), (28) 

= 6pi t i + (a - d)pi j2 - cp 2 ,2, (29) 

= 2(aCi,i - cCi, 2 ) 

+ 2^ 1 [ 7lil + (ZC-Zi2)Ci,i + (Zi2-ZC-l)77i i i], (30) 

= 2(6Ci, 2 -dC 2 , 2 ), (31) 

= &Ci,i + (a - d)Ci,2 - <2,2 

+ Kh^ + iZC-ZRK^ + iZR-ZC-l)^}, (32) 

= 2(0771,1 - c?7i i2 ) 

+ ( tv 2 - z - 1 ) [{ZR ~ zc ~ ^ " (33) 

= 2(H,2-dr ?2 , 2 ), (34) 

= 6771,1 + (a - d)?7i i2 - c?7 2j2 

+ ( n K -z- 1 ) " zc " 1)(Cl ' 2 " ^ 



(35) 



c = \ Wz^ ? ? ^ Qj ° ifc Q 7 ' ^ 

\ i j k c 

R = (j^EEE^^) , (37) 

where = (l//£!) and /n = #W(0) with #(0) = #( 2 )(0) = 0. 
2.3 Summary of our method 

The clustering coefficient C and the coupling connectivity R, which are given by Eqs. 
(36) and (37), respectively, play important roles in our DMA theory for SW networks. 



7 



The clustering coefficient C introduced in SW networks [3, 4], expresses a factor forming 
a cluster where the three sites i, j and k are mutually coupled. In contrast, the coupling 
connectivity R expresses a factor for a cluster where the two sites j and k are coupled 
to the third site i, but the sites j and k are not necessarily coupled [27]. It is noted in 
Eqs. (22)- (35) that there is no explicit dependence on the coupling randomness p, and it 
is only through parameters C and R in the mean-field equations. 

Figures 1(a) and 1(b) show the p dependences of C and R, respectively, for various 
Z values with N = 100. With increasing p from zero, C is decreased and approaches 
C = Z/N at p — 1. In contrast, R is monotonically increased with increasing p. 

Among the 12 correlations such as 7 KjA et al. given by Eqs. (15)-(18), nine correlations 
are independent because of the sum rule given by Eq. (20). In this study, we have 
chosen nine correlations of ^ K; \, Ck,a and p K ^\ as independent variables. Then the original 
2iV-dimensional stochastic DEs given by Eqs. (3) and (4) have been transformed to 11- 
dimensional deterministic DEs. Equations of motions for diffusive couplings given by Eqs. 
(22)-(35) are rather different from those for sigmoid couplings given by Eqs. (21)-(34) in 
[10], related discussion being given in Sec. 4. 

From a comparison of Eq. (24) with Eq. (27), we note that 

= for^/Z^O (38) 

= 71,1, for/V^^O (39) 

where Eq. (38) is nothing but the central-limit theorem describing the relation between 
fluctuations of local and average variables [Eqs. (15) and (19)]. In order to quantitatively 
discuss the synchronization, we first consider the quantity given by 

nt) = >^ E < MO - x iM 2 >= 2[7i,i(0 - PiM- (40) 

ij 

When all neurons are in the completely synchronous state, we get xu(i) = Xi(t) for all 
i, and then P(t) = in Eq. (40). On the contrary, in the asynchronous state, we get 
P(t) = 2(1 — l/N)^f 1:1 = Po(t) from Eq. (38). We have defined the synchronization ratio 
given by [10] 

^-iH ^-f" 1 )' 

which is and 1 for completely asynchronous (P = P ) and synchronous states (P = 0), 
respectively. 



8 



We define the time t max when S(t) takes its maximum value as 



t 



max 



{t | dS(t)/dt = 0,t in <t< t in + t w }. 



(42) 



The maximum value of S max [= S(t max )] depends on model parameters such as the cou- 
pling strength (K), the noise intensity {(3), the size of cluster (A), the coordination number 
(Z), and the coupling randomness (p), as will be discussed in the following section. 



3.1 Regular couplings 

We have adopted same parameters of 9 = 0.5, a = 0.5, t s = 10, A = 0.10, t in = 100 and 
T w = 10 as in [19], and the H function given by 



DMA calculations have been made by solving Eqs. (22)-(35) with the use of the fourth- 
order Runge-Kutta method with the time step of 0.01. We have performed also direct 
simulations by using the fourth-order Runge-Kutta method with the time step of 0.01. 
Results of direct simulations are averages of 1000 trials for Z < 20 and those of 100 trials 
otherwise noticed. All quantities are dimensionless. 

First we discuss the case of regular couplings (p = 0.0). Figs. 2(a), 2(b), 2(c), and 
2(d) show time courses of /ii, 7^, and respectively, with (3 = 0.005, K = 0.02, 
p = 0.0, N = 100 and Z = 10. Results of DMA expressed by solid curves are in good 
agreement with those of direct simulations depicted by dashed curves. Time courses of 
A*i) 7i,i Ci,i) an d Pi,i shown in Fig. 2(a)-2(d) are not so different from those for sigmoid 
couplings reported in Figs. 3(a)-3(d) of [10]. 

Fig. 3(a)-3(c) show time courses of S(t) calculated by DMA (solid curves) and direct 
simulations (dashed curves) for Z = 10, 50 and 99, whose magnitudes are increased with 
increasing Z. The maximum values of the synchronization ratio in the DMA are 0.0654, 
0.386 and 0.569 for Z = 10, 50 and 99, respectively, which shows a larger synchrony for 
larger Z. This is more clearly seen in Fig. 4(a) showing S max as a function of Z . Figure 
4(b) shows the Z dependences of 7^1, £1,1 and at t — t max with K = 0.02, (3 = 0.005 
and N = 100: filled and open marks express results of DMA and direct simulations, 
respectively. With increasing Z, 7^1 is significantly decreased while p± } i and are 
almost constant. This explains the larger synchrony Sf for larger Z, shown in Figs. 4(a). 



3 CALCULATED RESULTS 



H(xy — Xij) = Xij — Xij. 



(43) 
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The difference between the Z dependences of 7^1 and p^i is due to the fact that d'ji^/dt 
has a contribution from the second term of 2KZhi(£i t i — 7^1) in Eq. (24) while dpi^/dt 
has no such contributions in Eq. (27). Figure 4(b) shows that £1,1 also depends on Z 
because of the second term in Eq. (30). It is noted that because t max defined by Eq. 
(42) depends on Z in general, 7^1 at t — t max may show a weak Z dependence, as shown 
in Fig. 4(b) where t max = 107.16, 106.72, 106.46 and 105.96 for Z = 10, 20, 50 and 99, 
respectively. 

3.2 SW couplings 

Next we discuss the case of SW couplings by changing the coupling randomness p. Fig- 
ures 5(a), 5(b) and 5(c) show time courses of S(t) for p = 0.0, 0.1 and 1.0, respectively, 
calculated by the DMA (solid curves) and direct simulations (dashed curves). The maxi- 
mum values of the synchronization ratio S max in the DMA are 0.0654, 0.0694, and 0.0749, 
for p = 0.0, 0.1 and 1.0, respectively: S max is slightly increased with increasing p. This 
p dependence of S max is more clearly seen in Fig. 6(a) where S max is plotted against p 
for Z = 10. Figure 6(b) shows the p dependences of 7^1, £1,1 and at t — t max with 
K = 0.02, P = 0.005, iV = 100 and Z = 10: filled and open marks express results of 
DMA and direct simulations, respectively. With increasing p, 7^ is slightly decreased 
while pi 5 i is not changed. The origin of the difference between the p dependences of 71,1 
and pi t i is again due to the fact that dji^/dt has a contribution from the second term of 
2KZh 1 (( 1:1 — 7i 5 i) in Eq. (24) while dp lfl /dt has no such contributions in Eq. (27): the 
p dependence of 71,1 arises from (i ti which depends on p through network parameters of 
C and R in Eq. (30), as shown in Fig. 6(b). 

In Fig. 7, S max (p) normalized by its p = 0.0 value is plotted for various Z with 
N = 100, K = 0.02 and (3 = 0.005. Values of S max (p = 0.1)/ S max (p = 0.0) in the DMA 
are 1.061, 1.048, 1.0268 and 1.000 for Z = 10, 20, 30, and 50, respectively: an increase in 
S max is larger for smaller Z. 

4 CONCLUSION AND DISCUSSION 

Calculations in the preceding subsection show that when the coupling randomness p is 
introduced to regular networks, the synchronization may be slightly increased for diffusive 
couplings. This is in strong contrast with the result for sigmoid couplings in [10], which 
shows a decreased synchronization with increasing the coupling randomness. The main 
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origin of an increased synchronization for diffusive couplings may be their local negative 
feedback, as will be discussed in the followings. The diffusive coupling given by Eqs. (5) 
and (43) may be rewritten as 

l\ C) (t) = °ij ( x lj - x li) = K (J2 C V X ^ ~ k i X ^) ( 44 ) 

3 3 

where ki (= J2j is heterogeneous. We may show that the heterogeneity in the coor- 
dination number ki of Eq. (44) plays an important role in an increase of S max in SW 
networks. If we replace ki by its average of Z (=< ki >) in the feedback term of Eq. (44), 
it becomes 

I ( f\t)~K(£c ij x lj -Zx 1 ^. (45) 

3 

Filled and open circles in Fig. 8 denote S max calculated by using Eq. (44) with DMA and 
simulations, respectively, for N = 100, Z — 10, /3 — 0.005 and K = 0.02. Filled and open 
squares in Fig. 8 express S max calculated by using Eq. (45) with DMA and simulations, 
respectively, for the same parameters as mentioned above. Figure 8 clearly shows that 
heterogeneous negative-feedback term (—KkiXu) in Eq. (44) leads to a slightly increased 
synchronization whereas the homogeneous one (—KZxu) in Eq. (45) yields a decreased 
synchronization. 

Equation (44) may alternatively be rewritten as 

4 c) (t) = K ^2 dij xij, (46) 

3 

with 

dij Cij Sijki. (47) 

It is noted that the new adjacent matrix d^ given by Eq. (47) satisfies the relation given 
by 

£dy = 0. (48) 

3 

Nishikawa et al. [17] have studied the stability of synchronous states of coupled networks 
in which adjacent (Laplacian) matrix is assumed to satisfy the relation as given by Eq. 
(48). This implies that the coupling adopted in [17] is related to a diffusive process. 
From an analysis of the stability of synchronous state by the Lyapunov index, they have 
shown that the synchronization becomes more difficult in SW and SF networks with more 
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heterogeneity. Our calculation is expected not to be in contradict with theirs because they 
examine the criteria for the stability of synchronous oscillations while we have discussed 
the degree of the synchronization for an applied signal. It has been conventionally claimed 
that an increase in the synchronization arises from the shorten average distance L in SW 
networks [7, 11, 13, 14, 15]. However, equations of motions presented in Eqs. (22)-(35) 
(and those in [10]) do not include the term relevant to L of SW networks. 

There is also the difference between effects of heterogeneity for sigmoid and diffusive 
couplings. For sigmoid couplings [10], the effect of the heterogeneity of SW networks is 
included by a perturbation method with the term of 5cij [= Cij(p) — Qj(p = 0)] through 
new correlations functions of 0i and 2 [see Eq. (37) in [10]]. This has been made because 
the term of < SxuScij > appears in the process of calculating equations of motion, for 
example, of d^i/dt. In contrast, for the diffusive couplings, the counterpart term becomes 
< SxuSxij Scij >, which is in the higher order than < bx u bcij >. This shows that the 
effect of heterogeneity for diffusive couplings is weaker than that for sigmoid couplings: 
for the diffusive couplings its effect may be included by the p-dependent C and R in the 
mean-field approximation, while for sigmoid couplings it has to be taken into account by 
the perturbation method. The stronger heterogeneity for the sigmoid couplings yields a 
decrease in the synchronization when the heterogeneity is introduced. 

To summarize, we have discussed the synchronization in SW networks of spiking FN 
neurons with diffusive couplings, employing the semi- analytical DMA theory previously 
developed in [10]. A comparison of the results in this calculation with those for sigmoid 
couplings in [10] leads to the following conclusion. 

(1) When the average coordination number Z is increased, the synchronization S is in- 
creased both for sigmoid and diffusive couplings. We should note, however, that an 
increase in S is mainly made by an increase in for sigmoid couplings [Fig. 4(a) in 
[10]] while it is accomplished by a decrease in 7 1;1 for diffusive couplings [Fig. 4(b)]. 

(2) When the coupling randomness p is increased, the synchronization S is decreased by 
an introduced heterogeneity for sigmoid couplings, whereas for diffusive couplings, S may 
be slightly increased by their negative local feedback contribution which compensates its 
decrease caused by their heterogeneity. 

It is noted that an increase in the synchronization for diffusive couplings in the item (2) is 
not due to the shorten average distance in SW networks, against the conventional wisdom. 
The item (2) is consistent with the results in SW networks of the phase model with the 
coupling term of H(x — y) — sin(x — y), for which an increase in the synchronization with 
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increasing p has been reported [13, 14]. Items (1) and (2) imply that the synchronization 
of SW networks depends not only on the geometry of SW networks but also on details of 
couplings. In the present paper, we have neglected the transmission time delay. Because 
the average path length becomes shorter by added shortcuts [3, 4], the response speed is 
expected to be improved in SW networks with time delays. This is a great advantage of 
the SW networks though the synchronization may be not necessarily improved. Discus- 
sions in this paper and [10] have been confined to SW neural networks with symmetric 
(undirected) and unweighted couplings. Recently it has been shown that the synchroniza- 
tion in complex networks may be enhanced if their couplings are undirected and weighted 
[18]. It is interesting to apply our semianalytical approach to networks with directed and 
weighted couplings, which are realized in real complex networks. This subject is left as 
our future study. 
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APPENDIX A: Derivation of Eqs. (22)-(35) 

Substituting Eq. (14) to Eqs. (3)-(6), we get DEs for 5xu and 5x 2 i of a neuron i, 
given by (argument t is suppressed) 

dSxu 



dt 
d5x 



2j 



dt 



with 



dlPit) = K Y,Cij (h [8 Xlj (t) - 5 Xlj (t)} + h 3 [S Xlj (t) - Sx^it)} 3 ), (51) 



fi 8 Xli + h (5x 2 lt ~ 7i,i) + h 5x 3 lt - c 5x 2l + 5lt ] + (49) 
b 5x±j — d Sx 2 j, (50) 



where f e = (1/£\)F^ and h e = {l/l\)H^. DEs for the correlations are given by 



dt 
dt 



AT ^ 



\N 

(^EE% 



/ dSxxi \ , / dSx Ki \ 
ox Ki I — — I + I — — I 0X X i 



dt 



dt 



( d5x\j\ ( d5x 
Sx Ki I — 7-^ I + 



dt 



dt 



dp K ,\ I 1 \ - \ - 



/ c^aA , / d5x Kj 



dt 



(52) 
(53) 
(54) 



With the use of Eqs. (52)- (54), we may calculate DEs given by Eqs. (22)- (35). For 
example, terms including 51^ in d^i/dt, dCi^/dt and dp^i/dt become 



N 



E (SxuSPi 



2Khx 
N 



EE c u ((Sxii[Sx lj -Sx u \)) c , 



\ i 7 



2KZhi (Ci,i-7i,i), 



(55) 
(56) 
(57) 



^EE(W4 C) 



2KZh 1 
N 

2Kh x 



[71,1 + Z(C - i?)Ci,i + - - 1)771,1] 



EEE c ifc ((^ii[^ife-^u])) c , 



j j fe 



0, 



(58) 
(59) 

(60) 



In evaluating Eqs. (55)-(60), we have employed the relation given by Eq. (13): 



{{Sx Ki 5 x Xj )) c = 7 KjA Sij + C k ,a Cij + i] KyX (1 - Sij - Cij). 



(61) 
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Figure 1: The coupling randomness (p) dependence of (a) the clustering coefficient C and 
(b) the coupling connectivity R of SW networks for Z = 10, 20 and 50 with N = 100. 



Figure 2: (color online). Time courses of (a) /ii, (b) 7 1;1 , (c) and (d) pi ; i for (3 = 0.005, 
K = 0.02, iV = 100, Z = 10 and p = 0.0, solid and dashed curves denoting results 
of DMA and direct simulations, respectively. At the bottom of (a), an input signal is 
plotted. Vertical scales of (b), (c) and (d) are multiplied by factors of 10~ 4 , 10~ 5 and 
10~ 5 , respectively. 



Figure 3: (color online). Time courses of S(t) for (a) Z = 10, (b) Z = 50 and (c) Z = 99 
calculated by DMA (solid curves) and direct simulations (dashed curves) {(3 = 0.005, 
K = 0.02, N = 100 and p = 0.0). 



Figure 4: (color online). The average coordination number (Z) dependence of (a) S max 
and (b) 7^1 (circles), (triangles) and p^i (squares) at t — t max for (5 = 0.005, K = 0.02, 
iV = 100 and p = 0.0: filled and open marks denote results of DMA and direct simulations, 
respectively. 



Figure 5: (color online). Time courses of S(t) for (a) p = 0.0, (b) 0.1 and (c) 1.0 calculated 
by DMA (solid curves) and direct simulations (dashed curves) ((3 = 0.005, K = 0.02 and 
N = 100). 



Figure 6: (color online). The coupling randomness (p) dependence of (a) S max , and 
(b) 71,1 (circles), (triangles) and (squares) at t — t max for (3 = 0.005, K = 0.02, 
iV = 100 and Z = 10: filled and open marks denote results of DMA and direct simulations, 
respectively. 



Figure 7: (color online). The coupling randomness (p) dependence of S max (p) / S max (0) 
for Z = 10, 20 and 50 and with (3 = 0.005, K = 0.02, and N = 100. 



Figure 8: (color online). The coupling randomness (p) dependence of S max for the cou- 
plings KJ2j( c ij ~ $ijki) x ij (circles) and KJ2j( c ij ~ 8ijZ)x 1 j (squares) with N = 100, 
Z — 10, (3 — 0.005, and K = 0.02, filled and open marks denoting results of DMA and 
direct simulations, respectively (see text). 
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